#
# This code produces the experimental results in
# "You're Making Us Look Bad: Can Partisan Embarrassment
# Dampen Partisanship and Polarization?"
#

setwd("") # Note: set your own working directory here

#install.packages(c("tidyverse","stargazer","lmtest","AER")) # Uncomment if needed

library(tidyverse)
library(stargazer)
library(lmtest)
library(AER)

ex <- read.csv("study3.csv")

# Clean data
ex <- ex[-c(1:2),]
ex <- ex %>% filter(attn1 == "Bowling" & attn2 == "Orange")


# New variables
ex <- ex %>%
  mutate(democrat = as.numeric(pid == "Democrat" | 
                                 pid_lean == "Lean toward Democratic Party"),
         republican = as.numeric(pid == "Republican" |
                                   pid_lean == "Lean toward Republican Party"),
         pidstrength = case_when(
           pid_strength == "Not very strong ${q://QID2/ChoiceGroup/SelectedChoices}" ~ "Not very strong",
           pid_strength == "Strong ${q://QID2/ChoiceGroup/SelectedChoices}" ~ "Strong"),
         indy = as.numeric(pid == "Independent"),
         ideo7 = case_when(
           ideo == "Extremely liberal" ~ 1,
           ideo == "Liberal" ~ 2,
           ideo == "Slightly liberal" ~ 3,
           ideo == "Moderate" ~ 4,
           ideo == "Slightly conservative" ~ 5,
           ideo == "Conservative" ~ 6,
           ideo == "Extremely conservative" ~ 7),
         ba_plus = as.numeric(edu %in% c(
           "4-year college","Post-graduate or advanced degree")),
         nonwhite = as.numeric(race != "White"),
         treated = as.numeric(Condition == 2),
         control = as.numeric(Condition == 1))

# Dependent variables
ex <- ex %>%
  mutate(embarrassed = case_when(
    democrat == 1 & embarrassment_d == "Not at all embarrassed" ~ 0,
    democrat == 1 & embarrassment_d == "A little embarrassed" ~ 1,
    democrat == 1 & embarrassment_d == "Moderately embarrassed" ~ 2,
    democrat == 1 & embarrassment_d == "Very embarrassed" ~ 3,
    republican == 1 & embarrassment_r == "Not at all embarrassed" ~ 0,
    republican == 1 & embarrassment_r == "A little embarrassed" ~ 1,
    republican == 1 & embarrassment_r == "Moderately embarrassed" ~ 2,
    republican == 1 & embarrassment_r == "Very embarrassed" ~ 3),
    pid_importance = case_when(
      democrat == 1 & pid_important_d == "Not at all important" ~ 0,
      democrat == 1 & pid_important_d == "A little important" ~ 1,
      democrat == 1 & pid_important_d == "Moderately important" ~ 2,
      democrat == 1 & pid_important_d == "Very important" ~ 3,
      democrat == 1 & pid_important_d == "Extremely important" ~ 4,
      republican == 1 & pid_important_r == "Not at all important" ~ 0,
      republican == 1 & pid_important_r == "A little important" ~ 1,
      republican == 1 & pid_important_r == "Moderately important" ~ 2,
      republican == 1 & pid_important_r == "Very important" ~ 3,
      republican == 1 & pid_important_r == "Extremely important" ~ 4),
    social1 = case_when(
      democrat == 1 & social_d_1 == "Always" ~ 0,
      democrat == 1 & social_d_1 == "Most of the time" ~ 1,
      democrat == 1 & social_d_1 == "About half the time" ~ 2,
      democrat == 1 & social_d_1 == "Sometimes" ~ 3,
      democrat == 1 & social_d_1 == "Never" ~ 4,
      republican == 1 & social_r_1 == "Always" ~ 0,
      republican == 1 & social_r_1 == "Most of the time" ~ 1,
      republican == 1 & social_r_1 == "About half the time" ~ 2,
      republican == 1 & social_r_1 == "Sometimes" ~ 3,
      republican == 1 & social_r_1 == "Never" ~ 4),
    social2 = case_when(
      democrat == 1 & social_d_2 == "Always" ~ 0,
      democrat == 1 & social_d_2 == "Most of the time" ~ 1,
      democrat == 1 & social_d_2 == "About half the time" ~ 2,
      democrat == 1 & social_d_2 == "Sometimes" ~ 3,
      democrat == 1 & social_d_2 == "Never" ~ 4,
      republican == 1 & social_r_2 == "Always" ~ 0,
      republican == 1 & social_r_2 == "Most of the time" ~ 1,
      republican == 1 & social_r_2 == "About half the time" ~ 2,
      republican == 1 & social_r_2 == "Sometimes" ~ 3,
      republican == 1 & social_r_2 == "Never" ~ 4),
    social3 = case_when(
      democrat == 1 & social_d_3 == "Always" ~ 0,
      democrat == 1 & social_d_3 == "Most of the time" ~ 1,
      democrat == 1 & social_d_3 == "About half the time" ~ 2,
      democrat == 1 & social_d_3 == "Sometimes" ~ 3,
      democrat == 1 & social_d_3 == "Never" ~ 4,
      republican == 1 & social_r_3 == "Always" ~ 0,
      republican == 1 & social_r_3 == "Most of the time" ~ 1,
      republican == 1 & social_r_3 == "About half the time" ~ 2,
      republican == 1 & social_r_3 == "Sometimes" ~ 3,
      republican == 1 & social_r_3 == "Never" ~ 4),
    social4 = case_when(
      democrat == 1 & social_d_4 == "Always" ~ 0,
      democrat == 1 & social_d_4 == "Most of the time" ~ 1,
      democrat == 1 & social_d_4 == "About half the time" ~ 2,
      democrat == 1 & social_d_4 == "Sometimes" ~ 3,
      democrat == 1 & social_d_4 == "Never" ~ 4,
      republican == 1 & social_r_4 == "Always" ~ 0,
      republican == 1 & social_r_4 == "Most of the time" ~ 1,
      republican == 1 & social_r_4 == "About half the time" ~ 2,
      republican == 1 & social_r_4 == "Sometimes" ~ 3,
      republican == 1 & social_r_4 == "Never" ~ 4),
    shirt_public = case_when(
      democrat == 1 & support_d_1 == "Yes" ~ 1,
      republican == 1 & support_r_1 == "Yes" ~ 1,
      TRUE ~ 0),
    shirt_private = case_when(
      democrat == 1 & support_d_2 == "Yes" ~ 1,
      republican == 1 & support_r_2 == "Yes" ~ 1,
      TRUE ~ 0),
    donate_public = case_when(
      democrat == 1 & support_d_3 == "Yes" ~ 1,
      republican == 1 & support_r_3 == "Yes" ~ 1,
      TRUE ~ 0),
    donate_private = case_when(
      democrat == 1 & support_d_4 == "Yes" ~ 1,
      republican == 1 & support_r_4 == "Yes" ~ 1,
      TRUE ~ 0),
    volunteer_public = case_when(
      democrat == 1 & support_d_5 == "Yes" ~ 1,
      republican == 1 & support_r_5 == "Yes" ~ 1,
      TRUE ~ 0),
    volunteer_private = case_when(
      democrat == 1 & support_d_6 == "Yes" ~ 1,
      republican == 1 & support_r_6 == "Yes" ~ 1,
      TRUE ~ 0),
    pid_strength_post = case_when(
      pid_strength_today == "Not very strong ${q://QID24/ChoiceGroup/SelectedChoices}"~ 0,
      pid_strength_today == "Strong ${q://QID24/ChoiceGroup/SelectedChoices}" ~ 1)) 

ex$affpol_1 <- as.numeric(ex$affpol_1)
ex$affpol_2 <- as.numeric(ex$affpol_2)
ex$competence_1 <- as.numeric(ex$competence_1)
ex$competence_2 <- as.numeric(ex$competence_2)
ex$audience_1 <- as.numeric(ex$audience_1)
ex$audience_2 <- as.numeric(ex$audience_2)
ex$audience_3 <- as.numeric(ex$audience_3)

ex <- ex %>%
  rename(ft_post_gop = affpol_1,
         ft_post_dem = affpol_2,
         competence_gop = competence_1,
         competence_dems = competence_2) %>%
  mutate(ft_own_post = case_when(
    democrat == 1 ~ ft_post_dem,
    republican == 1 ~ ft_post_gop),
    ft_oppos_post = case_when(
      democrat == 1 ~ ft_post_gop,
      republican == 1 ~ ft_post_dem),
    competence_own = case_when(
      democrat == 1 ~ competence_dems,
      republican == 1 ~ competence_gop),
    inpartisans_inparty = case_when(
      democrat == 1 ~ other_evals_d_2,
      republican == 1 ~ other_evals_r_1),
    outpartisans_inparty = case_when(
      democrat == 1 ~ other_evals_d_1,
      republican == 1 ~ other_evals_r_2),
    indy_inparty = case_when(
      democrat == 1 ~ other_evals_d_3,
      republican == 1 ~ other_evals_r_3)) 

ex$outpartisans_inparty <- as.numeric(ex$outpartisans_inparty)
ex$inpartisans_inparty <- as.numeric(ex$inpartisans_inparty)
ex$indy_inparty <- as.numeric(ex$indy_inparty)

ex <- ex %>%
  mutate(ftownpost_01 = (ft_own_post - min(ft_own_post, na.rm=T))/(max(ft_own_post, na.rm=T)-min(ft_own_post, na.rm=T)),
         ftotherpost_01 = (ft_oppos_post - min(ft_oppos_post, na.rm=T))/(max(ft_oppos_post, na.rm=T)-min(ft_oppos_post, na.rm=T)),
         competence_01 = (competence_own - min(competence_own, na.rm=T))/(max(competence_own, na.rm=T)-min(competence_own, na.rm=T)),
         pid_importance01 = (pid_importance - min(pid_importance, na.rm=T))/(max(pid_importance, na.rm=T)-min(pid_importance, na.rm=T)),
         social1_01 = (social1 - min(social1, na.rm=T))/(max(social1, na.rm=T)-min(social1, na.rm=T)),
         social2_01 = (social2 - min(social2, na.rm=T))/(max(social2, na.rm=T)-min(social2, na.rm=T)),
         social3_01 = (social3 - min(social3, na.rm=T))/(max(social3, na.rm=T)-min(social3, na.rm=T)),
         social4_01 = (social4 - min(social4, na.rm=T))/(max(social4, na.rm=T)-min(social4, na.rm=T))
           )

#### Treatment effects ####
# Manipulation check
m1 <- lm(embarrassed ~ treated + democrat, data = ex)
bptest(m1)
m1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))

# ATE for own party FT
m2 <- lm(ftownpost_01 ~ treated + democrat, data = ex)
bptest(m2)

# Sub-group effects for own party FT
m2a <- lm(ft_own_post ~ treated, data = subset(ex, democrat == 1))
m2b <- lm(ft_own_post ~ treated, data = subset(ex, republican == 1))

# ATE for other party FT
m3 <- lm(ftotherpost_01 ~ treated + democrat, data = ex)
bptest(m3)
m3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC1"))

# Sub-group effects for other party FT
m3a <- lm(ft_oppos_post ~ treated, data = subset(ex, democrat == 1))
m3b <- lm(ft_oppos_post ~ treated, data = subset(ex, republican == 1))

# PID importance
m4 <- lm(pid_importance01 ~ treated + democrat, data = ex)
bptest(m4)

m4a <- lm(pid_importance01 ~ treated, data = subset(ex, democrat == 1))
m4b <- lm(pid_importance01 ~ treated, data = subset(ex, republican == 1))

# Competence
m5 <- lm(competence_01 ~ treated + democrat, data = ex)
bptest(m5)

m5a <- lm(competence_own ~ treated, data = subset(ex, democrat == 1))
m5b <- lm(competence_own ~ treated, data = subset(ex, republican == 1))

# Social polarization
m6 <- lm(social1_01 ~ treated + democrat, data = ex)
bptest(m6)

m6a <- lm(social1_01 ~ treated, data = subset(ex, democrat == 1))
m6b <- lm(social1_01 ~ treated, data = subset(ex, republican == 1))

m7 <- lm(social2_01 ~ treated + democrat, data = ex)
bptest(m7)

m7a <- lm(social2_01 ~ treated, data = subset(ex, democrat == 1))
m7b <- lm(social2_01 ~ treated, data = subset(ex, republican == 1))

m8 <- lm(social3_01 ~ treated + democrat, data = ex)
bptest(m8)

m8a <- lm(social3_01 ~ treated, data = subset(ex, democrat == 1))
m8b <- lm(social3_01 ~ treated, data = subset(ex, republican == 1))

m9 <- lm(social4_01 ~ treated + democrat, data = ex)
bptest(m9)

m9a <- lm(social4_01 ~ treated, data = subset(ex, democrat == 1))
m9b <- lm(social4_01 ~ treated, data = subset(ex, republican == 1))

# Campaign work
m10 <- lm(shirt_public ~ treated + democrat, data = ex)
bptest(m10)

m11 <- lm(shirt_private ~ treated + democrat, data = ex)
bptest(m11)

m12 <- lm(donate_public ~ treated + democrat, data = ex)
bptest(m12)
m12 <- coeftest(m12, vcov = vcovHC(m12, type = "HC1"))

m13 <- lm(donate_private ~ treated + democrat, data = ex)
bptest(m13)

m14 <- lm(volunteer_public ~ treated + democrat, data = ex)
bptest(m14)
m14 <- coeftest(m14, vcov = vcovHC(m14, type = "HC1"))

m15 <- lm(volunteer_private ~ treated + democrat, data = ex)
bptest(m15)

# PID strength post
m16 <- lm(pid_strength_post ~ treated + democrat, data = ex)

# Partisans' views of the in-party
# In-partisans' views of the in-party
in_in <- lm(inpartisans_inparty ~ treated + democrat, data = ex)
# Out-partisans' views of the in-party
out_in <- lm(outpartisans_inparty ~ treated + democrat, data = ex)
# Independents' views of the in-party
indy_in <- lm(indy_inparty ~ treated + democrat, data = ex)


#### Models interacted with PID strength ####
# ATE for own party FT
m2p <- lm(ftownpost_01 ~ treated*pidstrength + democrat, data = ex)
bptest(m2p)
m2p <- coeftest(m2p, vcov = vcovHC(m2p, type = "HC1"))
# ATE for other party FT
m3p <- lm(ftotherpost_01 ~ treated*pidstrength + democrat, data = ex)
bptest(m3p)
# PID importance
m4p <- lm(pid_importance01 ~ treated*pidstrength + democrat, data = ex)
bptest(m4p)
m4p <- coeftest(m4p, vcov = vcovHC(m4p, type = "HC1"))
# Competence
m5p <- lm(competence_01 ~ treated*pidstrength + democrat, data = ex)
bptest(m5p)
# Social polarization
m6p <- lm(social1_01 ~ treated*pidstrength + democrat, data = ex)
bptest(m6p)
m7p <- lm(social2_01 ~ treated*pidstrength + democrat, data = ex)
bptest(m7p)
m8p <- lm(social3_01 ~ treated*pidstrength + democrat, data = ex)
bptest(m8p)
m8p <- coeftest(m8p, vcov = vcovHC(m8p, type = "HC1"))
m9p <- lm(social4_01 ~ treated*pidstrength + democrat, data = ex)
bptest(m9p)
m9p <- coeftest(m9p, vcov = vcovHC(m9p, type = "HC1"))
# Campaign work
m10p <- lm(shirt_public ~ treated*pidstrength + democrat, data = ex)
bptest(m10p)
m10p <- coeftest(m10p, vcov = vcovHC(m10p, type = "HC1"))
m11p <- lm(shirt_private ~ treated*pidstrength + democrat, data = ex)
bptest(m11p)
m11p <- coeftest(m11p, vcov = vcovHC(m11p, type = "HC1"))
m12p <- lm(donate_public ~ treated*pidstrength + democrat, data = ex)
bptest(m12p)
m12p <- coeftest(m12p, vcov = vcovHC(m12p, type = "HC1"))
m13p <- lm(donate_private ~ treated*pidstrength + democrat, data = ex)
bptest(m13p)
m14p <- lm(volunteer_public ~ treated*pidstrength + democrat, data = ex)
bptest(m14p)
m14p <- coeftest(m14p, vcov = vcovHC(m14p, type = "HC1"))
m15p <- lm(volunteer_private ~ treated*pidstrength + democrat, data = ex)
bptest(m15p)
m15p <- coeftest(m15p, vcov = vcovHC(m15p, type = "HC1"))

#### IV analyses ####
# ATE for own party FT
m2b <- ivreg(ftownpost_01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m2b)
m2b <- coeftest(m2b, vcov = vcovHC(m2b, type = "HC1"))

# ATE for other party FT
m3b <- ivreg(ftotherpost_01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m3b)
m3b <- coeftest(m3b, vcov = vcovHC(m3b, type = "HC1"))

# PID importance
m4b <- ivreg(pid_importance01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m4b)

# Competence
m5b <- ivreg(competence_01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m5b)
m5b <- coeftest(m5b, vcov = vcovHC(m5b, type = "HC1"))

# Social polarization
m6b <- ivreg(social1_01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m6b)

m7b <- ivreg(social2_01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m7b)
m7b <- coeftest(m7b, vcov = vcovHC(m7b, type = "HC1"))

m8b <- ivreg(social3_01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m8b)

m9b <- ivreg(social4_01 ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m9b)

# Campaign work
m10b <- ivreg(shirt_public ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m10b)
m10b <- coeftest(m10b, vcov = vcovHC(m10b, type = "HC1"))

m11b <- ivreg(shirt_private ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m11b)

m12b <- ivreg(donate_public ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m12b)
m12b <- coeftest(m12b, vcov = vcovHC(m12b, type = "HC1"))

m13b <- ivreg(donate_private ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m13b)

m14b <- ivreg(volunteer_public ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m14b)
m14b <- coeftest(m14b, vcov = vcovHC(m14b, type = "HC1"))

m15b <- ivreg(volunteer_private ~ embarrassed + democrat | treated + democrat, data = ex)
bptest(m15b)
m15b <- coeftest(m15b, vcov = vcovHC(m15b, type = "HC1"))

#### Output ####
mods1 <- list(m2,m3,m4,m5,m6,m7,m8,m9) 
mods2 <- list(m10,m11,m12,m13,m14,m15)
mods3 <- list(m2b,m3b,m4b,m5b,m6b,m7b,m8b,m9b)
mods4 <- list(m10b,m11b,m12b,m13b,m14b,m15b)

stargazer(mods1, type = "text", out = "mods1.txt")
stargazer(mods2, type = "text", out = "mods1.txt")

stargazer(mods3, type = "text")
stargazer(mods4, type = "text")


#### Figures ####
#### Main coefficient plot ####
fx <- data.frame(
  Variable = c("FT Own Party","FT Other Party","Importance of PID",
               "Party is competent","Do favors","Watch property",
               "Ask personal things","Talk about politics",
               "Wear shirt (public)","Wear shirt (private)",
               "Donate (public)","Donate (private)",
               "Volunteer (public)","Volunteer (private)"),
  Effect = c(m2$coefficients[2], m3[,1][2], m4$coefficients[2], m5$coefficients[2],
             m6$coefficients[2], m7$coefficients[2], m8$coefficients[2], m9$coefficients[2],
             m10$coefficients[2], m11$coefficients[2], m12[,1][2], m13$coefficients[2],
             m14[,1][2], m15$coefficients[2]),
  SE = c(coef(summary(m2))[, "Std. Error"][2], m3[,2][2],
         coef(summary(m4))[, "Std. Error"][2], coef(summary(m5))[, "Std. Error"][2],
         coef(summary(m6))[, "Std. Error"][2], coef(summary(m7))[, "Std. Error"][2],
         coef(summary(m8))[, "Std. Error"][2], coef(summary(m9))[, "Std. Error"][2],
         coef(summary(m10))[, "Std. Error"][2], coef(summary(m11))[, "Std. Error"][2],
         m12[,2][2], coef(summary(m13))[, "Std. Error"][2],
         m14[,2][2], coef(summary(m15))[, "Std. Error"][2])
)

interval95 <- -qnorm((1-0.95)/2)

fx$Variable <- factor(fx$Variable, levels = fx$Variable[order(fx$Effect)])

coefs <- ggplot(fx)
coefs <- coefs + geom_hline(yintercept = 0, color = gray(1/2), lty = 2)
coefs <- coefs + geom_linerange(aes(x = Variable, ymin = Effect - SE*interval95, ymax = Effect + SE*interval95), 
                                lwd = 1, position = position_dodge(width = 1/2))
coefs <- coefs + geom_pointrange(aes(x = Variable, y = Effect, ymin = Effect - SE*interval95, ymax = Effect + SE*interval95),
                                 lwd = 1/2, position = position_dodge(width = 1/2), shape = 21, fill = "BLACK")
coefs <- coefs + coord_flip() + theme_classic()
coefs <- coefs + ylab("Treatment Effect") + xlab("Dependent Variable")
ggsave("Figure 3.png", dpi = 600)

#### PID Strength coefficient plot ####
fx <- data.frame(
  Variable = c("FT Own Party","FT Other Party","Importance of PID",
               "Party is competent","Do favors","Watch property",
               "Ask personal things","Talk about politics",
               "Wear shirt (public)","Wear shirt (private)",
               "Donate (public)","Donate (private)",
               "Volunteer (public)","Volunteer (private)"),
  Effect = c(m2p[,1][5], m3p$coefficients[5], m4p[,1][5], m5p$coefficients[5],
             m6p$coefficients[5], m7p[,1][5], m8p[,1][5], m9p[,1][5],
             m10p[,1][5], m11p[,1][5], m12p[,1][5], m13p$coefficients[5],
             m14p[,1][5], m15p[,1][5]),
  SE = c(m2p[,2][5], coef(summary(m3p))[, "Std. Error"][5], 
         m4p[,2][5], coef(summary(m5p))[, "Std. Error"][5],
         coef(summary(m6p))[, "Std. Error"][5], m7p[,2][5],
         m8p[,2][5], m9p[,2][5],
         m10p[,2][5], m11p[,2][5],
         m12p[,2][5], coef(summary(m13p))[, "Std. Error"][5],
         m14p[,2][5], m15p[,2][5])
)

interval95 <- -qnorm((1-0.95)/2)


fx$Variable <- factor(fx$Variable, levels = fx$Variable[order(fx$Effect)])

coefs <- ggplot(fx)
coefs <- coefs + geom_hline(yintercept = 0, color = gray(1/2), lty = 2)
coefs <- coefs + geom_linerange(aes(x = Variable, ymin = Effect - SE*interval95, ymax = Effect + SE*interval95), 
                                lwd = 1, position = position_dodge(width = 1/2))
coefs <- coefs + geom_pointrange(aes(x = Variable, y = Effect, ymin = Effect - SE*interval95, ymax = Effect + SE*interval95),
                                 lwd = 1/2, position = position_dodge(width = 1/2), shape = 21, fill = "BLACK")
coefs <- coefs + coord_flip() + theme_classic()
coefs <- coefs + ylab("Treatment Effect") + xlab("Dependent Variable")

#### IV coefficient plot ####
fx <- data.frame(
  Variable = c("FT Own Party","FT Other Party","Importance of PID",
               "Party is competent","Do favors","Watch property",
               "Ask personal things","Talk about politics",
               "Wear shirt (public)","Wear shirt (private)",
               "Donate (public)","Donate (private)",
               "Volunteer (public)","Volunteer (private)"),
  Effect = c(m2b[,1][2], m3b[,1][2], m4b$coefficients[2], m5b[,1][2],
             m6b$coefficients[2], m7b[,1][2], m8b$coefficients[2], m9b$coefficients[2],
             m10b[,1][2], m11b$coefficients[2], m12b[,1][2], m13b$coefficients[2],
             m14b[,1][2], m15b[,1][2]),
  SE = c(m2b[,2][2], m3b[,2][2], coef(summary(m4b))[, "Std. Error"][2], m5b[,2][2],
         coef(summary(m6b))[, "Std. Error"][2], m7b[,2][2], coef(summary(m8b))[, "Std. Error"][2], 
         coef(summary(m9b))[, "Std. Error"][2], m10b[,2][2], coef(summary(m11b))[, "Std. Error"][2],
         m12b[,2][2], coef(summary(m13b))[, "Std. Error"][2], m14b[,2][2], m15b[,2][2])
)

interval95 <- -qnorm((1-0.95)/2)

fx$Variable <- factor(fx$Variable, levels = fx$Variable[order(fx$Effect)])

coefs <- ggplot(fx)
coefs <- coefs + geom_hline(yintercept = 0, color = gray(1/2), lty = 2)
coefs <- coefs + geom_linerange(aes(x = Variable, ymin = Effect - SE*interval95, ymax = Effect + SE*interval95), 
                                lwd = 1, position = position_dodge(width = 1/2))
coefs <- coefs + geom_pointrange(aes(x = Variable, y = Effect, ymin = Effect - SE*interval95, ymax = Effect + SE*interval95),
                                 lwd = 1/2, position = position_dodge(width = 1/2), shape = 21, fill = "BLACK")
coefs <- coefs + coord_flip() + theme_classic()
coefs <- coefs + ylab("Induced Embarrassment") + xlab("Dependent Variable")


